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Abstract:  It  is  known  that  scattering  by  particulates  within  natural  waters  is 
the  main  cause  of  the  blur  in  underwater  images.  Underwater  images  can  be 
better  restored  or  enhanced  with  knowledge  of  the  point  spread  function 
(PSF)  of  the  water.  This  will  extend  the  performance  range  as  well  as  the 
information  retrieval  from  underwater  electro-optical  systems,  which  is 
critical  in  many  civilian  and  military  applications,  including  target  and 
especially  mine  detection,  search  and  rescue,  and  diver  visibility.  A  better 
understanding  of  the  physical  process  involved  also  helps  to  predict  system 
performance  and  simulate  it  accurately  on  demand.  The  presented  effort 
first  reviews  several  PSF  models,  including  the  introduction  of  a  semi- 
analytical  PSF  given  optical  properties  of  the  medium,  including  scattering 
albedo,  mean  scattering  angles  and  the  optical  range.  The  models  under 
comparison  include  the  empirical  model  of  Duntley,  a  modified  PSF  model 
by  Dolin  el  aly  as  well  as  the  numerical  integration  of  analytical  forms  from 
Wells,  as  a  benchmark  of  theoretical  results.  For  experimental  results,  in 
addition  to  that  of  Duntley,  we  validate  the  above  models  with  measured 
point  spread  functions  by  applying  field  measured  scattering  properties  with 
Monte  Carlo  simulations.  Results  from  these  comparisons  suggest  it  is 
sufficient  but  necessary  to  have  the  three  parameters  listed  above  to  model 
PSFs.  The  simplified  approach  introduced  also  provides  adequate  accuracy 
and  flexibility  for  imaging  applications,  as  shown  by  examples  of  restored 
underwater  images. 
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1.  Introduction 

Image  formation  is  usually  affected  by  both  the  optical  system  design  including  lenses  and 
recording  devices,  as  well  as  transfer  processes  associated  within  the  formation  medium. 
Applications  to  remove  the  effects  of  the  medium  have  largely  been  done  in  atmospheric 
situations,  where  astronomical  and  reconnaissance  needs  demand  a  better  understanding  and 
solution.  This  typically  involves  the  effects  due  to  atmospheric  turbulence.  In  natural  waters, 
although  the  turbulence  does  impose  its  limitations  under  certain  conditions,  especially  at 
longer  range[l,  2],  the  largest  contribution  to  degraded  images  comes  from  suspended 
particulates.  Their  size  and  index  of  refraction  changes  cause  blurring  in  the  resulting 
imagery,  thus  limiting  visibility  range. 

In  general,  scattering  properties  of  the  medium  determine  the  outcome  of  the  image 
transmission.  In  ocean  and  lake  environment  optics,  the  scattering  properties  are  conveniently 
described  and  measured  by  the  scattering  coefficient  ( b ),  which  determines  the  probability  that 
a  photon  will  be  scattered  away  from  its  original  traveling  direction  per  unit  length  by  the 
medium  molecules,  constituents  within  the  medium  (ie  particles  [3]),  and  turbulence  [2].  The 
scattering  parameter  (6)  is  an  integration  of  the  volume  scattering  or  phase  function,  /?,  which 
details  such  probabilities  by  the  relative  directions  of  incoming  and  out-going  photons[3]. 
These  inherent  optical  properties  (IOP),  although  measured  frequently  due  to  their  important 
applications  in  ocean  optics,  especially  in  remote  sensing,  have  not  been  applied  to 
underwater  imaging  issues  directly,  in  part  because  they  inherently  only  reflect  the  effect  of 
the  single  scattering. 

A  more  intuitive  and  direct  measure  in  imaging  is  the  point  spread  function  (PSF),  which 
gives  the  system  response  to  a  point  source,  and  thus  includes  the  effect  of  multiple  scattering 
[3,  4].  It  is  the  ideal  parameter  to  study  image  transmission,  optical  sounding  [5],  and  retrieval 
of  optical  properties  [6].  Simply  put,  an  image  of  an  object  is  the  combination  of  the  original 
signal,  f(x,y)y  convolved  with  the  entire  system  point  spread  function  h(x,y)9  adding  noise 
n(x,y), 

g(x,y)  =  f(x,y)®h(x,y)  +  n(x,y),  (1) 

The  system  response  includes  those  from  both  the  imaging  system  itself,  as  well  as  the 
effects  of  the  medium  (water  in  our  case).  With  known  characteristics  of  the  imaging  system 
and  correct  modeling  of  the  medium,  theoretically  it  is  possible  to  fully  recover  the  original 
signal  by  reversion  or  deconvolution[7].  Mathematically  PSF  is  equivalent  to  the  beam  spread 
function  (BSF)  [8]  which  can  be  modeled  and  measured  more  easily.  From  (1),  one  can  see 
that  with  the  knowledge  of  g  and  /?,  the  original  image  / can  be  restored,  especially  when  the 
noise  can  be  reduced  with  improved  system  setup  or  optimization^].  Additionally,  with  the 
knowledge  of  the  PSF,  one  can  easily  simulate  the  system  outcome.  This  is  valuable  to  system 
design,  performance  prediction,  as  well  as  underwater  scene  simulations,  as  it  is  based  on  the 
accuracy  of  the  physical  model  of  the  system  (medium  included). 
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2.  Selected  previous  PSF  models 

Although  established  almost  40  years  ago,  the  most  thorough  empirical  point  spread  model  to 
date  is  still  that  of  Duntley  [10].  Duntley  reported  extensive  lab  measurements  of  BSFs  of 
simulated  ocean  waters  with  remarkably  different  optical  lengths  (0.5  to  21),  and  summarized 
his  Findings  in  a  single,  albeit  complex  empirical  relationship: 

H(0)  =  \OAC0B  (2) 

P  2ks\x\0 


where  parameters  £  =  1-2(10"°),  =  1.260 -0.375r° 7,0+0 489/f -(1.378 +  0.053^)1 0~(0 268+0 088)r, 


£T  =  (13.75- 0.50 If) - (0.626 - 0.357f  )r  +  (0.01 258  +  0.00354f )r2 ,  C7  =  — [1— (—  +  £)2'}],  and 

3  E 

D  =  (0.01 8 +  0.01 1£  + 0.001 725r)r .  Note  in  the  above  equations  that  the  optical  length  is 
defined  as  r  =  cr  (c  being  the  total  attenuation  coefficient)  while  the  single  scattering  albedo 
is  given  as  ^  =  b/c .  g  =  1/(1  -  o^)  is  used  to  simplify  the  form.  H  is  the  spectral  irradiance 


and  is  normalized  by  the  laser  power  P,  which  gives  the  BSF.  To  add  to  the  complexity,  one 
should  realize  that  the  above  regression  results  are  in  degrees,  rather  than  radians.  It  is  of  little 
surprise  that  the  original  publication  had  to  resolve  to  errata  to  get  the  exact  form  correct, 
while  another  citation  also  contained  misinformation  likely  due  to  printing  errors  [11]. 

A  different  approach  was  used  by  Dolin  et  al  [12],  where  numerical  approximations  were 
applied  in  his  analytical  formulation  of  the  PSF.  We  simplify  this  model  by  assuming  little 
backscattering  compared  to  the  total  scattering.  Briefly,  a  modified  BSF  E(0,r)  can  be 
expressed  via  PSF  kernel  G,  in  the  order  of  direct  beam,  single  and  multiple  scattering  terms 
as: 


E(0yr)=  P0(q  /  r  )2  G(0q1Th)  exp  (~aef  r  ), 


(3) 


G{Oq,Th)=  exp(-r„ ) + 0.525^- exp(-2.60,° 7  -  rh )  +  &-[2 - (1  +  Tb ) exp(-rk )] 

7tu„  u ,  LK 

</  ? 


(4) 


6.857  - 1 .5737r,  +  0. 143rt2  -  6.027  10  3  •  rj,  + 1 .069  1  (T*  rt4 

1-0.1 869r„  + 1 .97  1 0-2  r]  - 1 .3 1 7  •  1 0'3  •  rt3  +  4.3 1  1 0'5  r*  ’ 


B  _  0.469 -7.41- 10'2r,  +2.78  10  3rt2 +9.6- IQ'5  rt3 

Pl~  1-9.16  10-2r» -6.07  10'3  r2  +8.33  10^  r2 

6.27 -0.723r,,  + 5.82  10~2-f2 
~  1  -  0.072rft  +  6.3  •  1 0‘J  •  r2  +  9.4  1 0"  •  r\  ' 

0q=q9  ,Th=br  =  a)T, 


(6) 

(7) 

(8) 


27t^G(dq,Th)0qddq  =  1, 


]s{6q)0qdeq  =  \- 


1- 


cos£ 


1-2  b, 


b  _ 


-4  1- 


cos£ 


1-2  bh 


bh=bh/b. 


(9) 

(10) 


Notice  that  a  dimensionless  scaling  factor  q  is  used,  which  can  be  shown  to  approximate 
the  mean  scattering  angle  [12].  P0  and  S  relate  to  the  direct  beam  initial  power  and  can  be 
omitted  when  normalizations  are  applied  to  scattered  photons.  The  above  simplified  result  is 
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under  the  assumption  that  the  absorption  coefficient  a»bb ,  the  backscattering  coefficient,  and 
b»bb ,  such  that  b±2bb~b ,  and  a±2bb~a.  This  holds  true  in  most  conditions  applied  and 
introduces  little  variation  when  tested.  This  is  not  surprising  as  our  focus  here  is  on  small 
forward  scattering  angles.  Once  again,  the  complexity  of  this  model  is  evident  and  even 
exceeds  that  of  Duntley’s.  However,  it  does  offer  a  better  approximation,  as  we  shall  examine 
later. 

Lastly,  using  field  collected  data,  Voss  put  forth  a  simple-looking  empirical  PSF 
expression  [11],  that  fits  3  different  types  of  oceanic  waters  (TOTO,  Pacific  and  Sargasso 
Sea).  The  PSF  covered  the  range  between  4  and  100  mrad  and  was  given  by: 

PSF(O)  =  Bfim ,  (11) 

where  -m  is  the  slope  of  log(PSF)  -  log(0),  and  B/  is  a  constant.  Here  m  is  not  a  constant,  but 
rather  a  function  of  lOP  and  r.  The  above  formula  fits  all  data  to  less  than  15%  error.  For 
imaging  restoration  applications,  especially  considering  optimization  in  convolution[9],  this  is 
an  acceptable  start.  The  m  values  ranged  from  0.4  to  2.0  as  a  function  of  r(0  to  10).  Other 
than  graphical  results,  there  is  no  definitive  relationship  provided  to  allow  comparison  of  this 
method  to  other  approaches.  Nonetheless,  this  is  a  very  encouraging  result,  since  such  simple 
relationships  can  be  rather  beneficial  to  imaging  needs,  especially  when  high  frame  rate  per 
pixel  calculations  are  needed.  The  aim  of  this  paper  is  to  explore  this  form  further  by 
examining  the  PSF  relationship  to  commonly  measured  IOPs  in  search  of  an  even  simpler 
solution  for  real-time  imaging  needs. 

3.  A  simplified  PSF  model 

A  simple  form  of  the  PSF  is  highly  desirable,  especially  when  routinely  measured  optical 
properties  can  be  used  as  input.  We  introduced  such  a  model  previously!  1 3],  based  on  the 
analytical  form  by  Gordon[14].  For  brevity  only  key  steps  are  outlined  here.  We  first  note  that 
the  PSF  can  be  approximated  by  examining  the  amount  of  non-disturbed  light  variation  within 
a  cone  of  half  angle  0  at  range  r  [14,  15],  therefore 

PSF(6,r)  =  —  -  (12) 


l  V/t 

<t>(0,r)  =  exp[-(c- J  J pd(adt)r]' 


PSF(e)  =  e’r\fi( 

A  *  * 


(13) 

(14) 


Direct  integration  of  (14)  can  be  carried  out  for  a  known  analytical  volume  scattering 
function.  The  parameter  t  is  used  for  the  convenience  of  integration  and  can  be  interpreted  as 
the  portion  of  the  scattering  angle  up  to  6.  There  are  many  possible  and  defendable  choices  of 
p  of  natural  waters.  Here  we  take  a  simple  formulation  from  Wells  [  1 5],  by  using  a  phase 
function  in  the  following  form: 


bOn 


2n(pl+e2y 


15) 


where  60  relates  to  the  mean  scattering  angle.  Substituting  (15)  into  (14)  for  integration  will 
give  inverse  power  relationships  to  the  scattering  angle  6.  With  the  help  of  numerical 
integrations,  we  arrive  at  the  following  empirical  form: 


PSF(0)  =  K(60 )^— 
0  2  nff" 


=  *«U 


2nd" 


(16) 
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where  K  is  a  constant  and  does  not  affect  imaging  needs  since  relative  units  of  PSF  are  used, 
and  m=l/(Oo-2T0o.  The  later  term  is  necessary  to  include  higher  scattering  orders.  This  is 
inline  with  numerical  calculations  over  different  parameter  ranges.  Comparisons  to  empirical 
and  measured  formations  can  now  be  made. 

4.  Comparison  and  validation  of  PSF  models 

We  proceed  by  comparing  the  following  PSF  models:  1)  our  simplified  model;  2)  Duntley’s 
empirical  relationship  from  simulated  seawaters;  3)  a  modified  Dolin’s  PSF  model  from 
analytical  results;  4)  numerically  calculated  PSF  using  Wells  analytical  method;  and  5)  a  PSF 
derived  from  field  scattering  measurements  and  applying  Monte  Carlo  (MC)  methods.  We  see 
that  the  second  and  third  methods  are  straight-forward.  The  fourth  approach,  found  in  other 
publications  [16,  17],  offers  a  means  to  use  ground-truth,  to  accurately  simulate  PSFs  from 
optical  scattering  properties.  The  details  of  this  approach  can  be  found  in  a  related  paper  [4], 
but  a  brief  description  is  offered  below. 

This  method  uses  field  measured  lOPs  including  volume  scattering  functions  to 
numerically  calculate  BSF,  thus  PSF,  using  MC.  As  a  proven  technique,  MC  is  essentially  a 
numerical  experiment  of  the  physical  process  under  consideration,  thus  providing  accurate 
results  although  it  does  require  extensive  computation  as  the  results  are  statistical  in  nature 
and  each  set  of  parameters  must  be  run  separately.  In  fact,  the  MC  method  has  been  used 
previously  to  test  the  range  dependence  of  the  Wells  theory[18].  Given  the  IOPs  (phase 
function,  absorption  and  scattering  coefficients),  the  radius  of  the  sphere  R,  and  an  angular 
resolution  AO,  the  code  produces  the  irradiance  for  each  angular  bin  over  the  entire  sphere. 

Data  were  obtained  during  an  April-May  2006  NATO  trial  experiment  in  Panama  City, 
Florida.  The  amount  of  scattering  and  absorption  were  controlled  by  introducing  Maalox  and 
absorption  dye  respectively  into  a  tap-water  filled  pool.  In-water  optical  properties  during  the 
experiment  were  measured  and  included  the  absorption  and  attenuation  coefficients  (Wetlabs 
ac-9),  particle  size  distributions  (Sequoia  Scientific  LISST-100),  and  volume  scattering 
functions  (multi-spectral  volume  scattering  meter  or  MVSM). 


Tabic  1 .  Optical  parameters  used  to  compare  different  PSF  models  as  shown  in  Figs.  1-2. 


Date 

optical 

measured 

scattering 

length  (x) 

c  (m ') 

albedo  (coo) 

4/27/2006 

3.0 

0.78 

0.87 

4/27/2006 

5.5 

0.78 

0.87 

5/1/2006 

3.7 

1.38 

0.89 

5/1/2006 

7.3 

1.38 

0.89 

5/2/2006 

2.5 

0.91 

0.64 

5/2/2006 

5.5 

0.91 

0.64 

A  comparison  of  PSFs  from  different  methods  outlined,  namely,  the  current  study,  MC 
simulated  PSF  from  measured  MVSM  phase  functions.  Wells’  numerical  results,  Dolin’s 
model,  and  Duntley’s  are  shown  in  Figs.  1  and  2,  and  parameters  listed  in  Table  1.  These  are 
normalized  to  lOOmrad  values.  The  current  approach  agrees  well  with  measurements  as  well 
as  empirical  models  under  the  different  testing  conditions  shown.  There  are  only  small 
differences  except  at  larger  scattering  angles  with  the  MVSM,  Wells  and  Dolin.  Different 
optical  range  and  scattering  conditions  are  compared  in  Fig.  1,  under  similar  scattering 
characteristics  (00=0.03),  which  would  correspond  to  <7=7.6  in  Dolin’s  model.  We  can  see  that 
Dolin  and  Wells  models  match  those  of  measured  MVSM  very  well  under  all  conditions, 
which  is  not  surprising  as  they  applied  few  approximations  during  derivation. 
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Fig.  1.  Comparison  of  PSFs  from  different  methods  from  a-f:  “Current”  (current  study), 
“Duntley”  and  “Dolin”  (empirical  relationship  in  text),  and  “MVSM”  (Monte  Carlo  simulated 
PSF  based  on  measured  volume  scattering  function  during  NATO  experiment)  and  Wells 
numerical  integrated  PSF.  Normalized  to  100  mrad.  See  text  for  details. 


At  larger  angles,  there  are  visible  differences  between  the  current  model  results  and  those 
of  Wells  and  Dolin,  however,  it  still  approximates  Duntley’s  empirical  results  nicely.  These 
differences  only  have  small  effects  on  image  restoration[6],  because  of  the  relatively  small 
effects  (log  scale)  in  grayscale  pixel  values  (ie,  1  or  2  out  of  255,  occurring  several  pixels 
away  from  the  center  of  PSF).  Initial  restoration  testing  with  images  obtained  during  the 
NATO  experiment  show  little  differences  using  any  of  these  PSFs  under  typical  conditions 
tested[17].  If  simplicity  or  computation  time  is  of  less  concern,  the  differences  with  the 
current  model  at  larger  angles  can  be  effectively  reduced  by  applying  a  “damping  factor”, 
such  as  exp(-50‘ 2)  to  the  right  side  of  Eq.  (16),  or  apply  Wells’  or  Dolin’s  model  directly.  One 
may  notice  from  the  figure  that  PSFs  from  MVSM  are  almost  always  higher  than  Duntley’s. 
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This  is  likely  due  to  the  differences  in  the  scattering  agents  involved,  as  the  particle  sizes  will 
be  different  because  of  the  different  preparation  process  (recall  that  Duntley  used  high  purity 
filtered  water  as  their  experiment  base  in  the  tank,  before  adding  similar  absorption  and 
scattering  agents).  In  fact,  the  differences  in  PSFs  can  be  as  high  as  an  order  of  magnitude 
under  certain  conditions  at  small  angles  (eg.  0*7=0.87,  T=5.5),  or  ±1  in  PSF  log  slope. 
Examination  of  the  parameter  related  to  Go  suggests  that  a  higher  value  (0.06)  is  more 
appropriate  in  order  to  match  Duntley’s  PSF,  corresponding  to  less  forward  scattering 
contribution  in  Duntley’s  case  (Fig.  2).  This  is  not  surprising  since  despite  the  complexity  of 
Duntley’s  model,  it  depends  on  only  two  parameters,  namely  the  optical  range  and  scattering 
albedo.  It  lacks  the  ability  to  alter  the  shape  of  the  volume  scattering  function.  Variations  by 
Go  (Fig.  2)  demonstrated  the  flexibility  of  current  model,  and  also  suggest  that  the  reason 
behind  Voss’  choice  of  tyf=0.06  [11]  may  likely  be  the  result  of  different  particle  types  thus 
the  shape  of  scattering  functions  involved.  Also  note  that  the  PSF  from  the  Wells  model  is 
more  responsive  than  that  of  Dolin.  Also,  it  is  worth  pointing  out,  that  the  current  approach 
matches  those  from  Duntley’s  well  to  at  least  r=15,  but  with  much  less  computation  required. 


Fig.  2.  An  example  of  FSF  comparison  under  different  0O  values.  Notice  that  the  variation  of  0O 
docs  not  affect  Duntley  or  MVSM  and  all  rest  parameters  of  other  models  arc  the  same  as  in 
Fig.l.  Normalized  to  100  mrad. 


The  above  derived  PSF  model  (Eq.  16)  was  further  validated  when  applied  to  restore 
blurred  underwater  images.  A  Secchi  disk  image  taken  by  Navy  diver’s  Nikon  Coolpix  digital 
camera  at  5  ft  away  during  an  exercise  to  test  underwater  diver  visibility  models  (August  21, 
2002,  Mississippi  Bight,  c=1.5m'1)  was  used.  When  the  correct  PSF  is  applied,  the  original 
image  showed  significant  improvements,  both  in  details  and  contrast  [Fig.  3(b)].  When  the 
same  formula  is  used  for  the  PSF  but  assuming  a  longer  optical  length  (x=10  instead  of  2.3), 
we  notice  actual  degradation  in  details  due  to  over-correction,  although  contrast  was  still 
improved.  Both  restorations  were  done  using  an  iterative  Lucy-Richardson  method[7].  When 
a  blind  deconvolution  approach  [7]  was  tested  with  an  intentionally  different  PSF  (Gaussian  in 
this  case),  the  result  was  apparently  worse,  even  though  traditionally  (ie  in  astronomy 
applications)  such  algorithms  are  not  sensitive  to  the  shape  of  PSFs  [7,  9].  These  examples 
farther  demonstrate  the  importance  and  validity  of  the  correct  PSF  model  to  underwater  image 
restorations. 
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Fig.  3.  A  diver  captured  image  and  restorations  with  different  PSFs.  Image  (top  left,  3a)  was 
taken  during  August  21,  2002  diving  at  Mississippi  Bight  with  total  attenuation  c=1.5m  1  at  5ft 
distance.  3b  (top  right)  shows  the  restored  image  using  the  correct  PSF  model.  3c  (bottom  left) 
shows  the  image  restored  assuming  a  longer  distance.  3d  (bottom  right)  shows  the  result  of  a 
blind  deconvolution  with  a  Gaussian  PSF. 


5.  Conclusion 

Several  types  of  PSF  models  of  natural  waters  including  empirical,  analytical  and 
measurement  results  under  different  conditions  are  compared,  along  with  the  introduction  of  a 
simplified  semi-analytical  form  which  shows  a  close  approximation  to  the  above  methods. 
Measured  PSFs  were  obtained  using  MVSM-measured  phase  functions  and  MC  simulations. 
The  results  show  the  PSFs  from  a  modified  Dolin  and  Wells  are  very  accurate  even  at  larger 
angles  and  longer  optical  ranges.  It  is  shown  that  three  parameters  are  needed  to  adequately 
model  the  system  response  function  of  the  natural  waters.  This  explains  in  part  the 
inflexibility  of  Duntley’s  empirical  relationship.  Although  small  differences  exist  between  the 
simplified  model  and  other  approaches,  they  are  not  significant  as  far  as  image  restoration  is 
concerned.  The  simpler  form  of  the  model  can  be  beneficial  to  various  imaging  needs  in 
natural  environments,  especially  to  real-time  processing  and  scene  generation  applications. 
The  fact  that  different  scattering  agents  can  result  in  an  order  of  magnitude  difference  in  PSFs 
between  measured  MC-MVSM  and  Duntley’s  at  small  angles,  which  affects  the  slope  of  the 
PSF  significantly,  suggests  additional  studies  are  needed  to  better  determine  the  effects  of 
particle  types  on  the  shape  of  scattering  functions. 
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